pool_mfx_summary <- function(list_mfx_summaries) {
  ## Sub-functions
  pool_imputations <- function(list_mfx_summaries) {
    est <- sapply(list_mfx_summaries, function(x)
      x$AME)
    se <- sapply(list_mfx_summaries, function(x)
      x$SE)
    
    if (!is.matrix(est)) {
      est <- t(est)
      se <- t(se)
    }
    
    ## Point estimates
    est_pooled <-  apply(est, 1, mean)
    
    ## Standard errors
    pool_ses <- function(est, se) {
      m <- ncol(se)
      V_W <- apply(se, 1, function(x)
        mean(x ^ 2))
      V_B <- (1 - m) * apply(est, 1, function (x)
        sum(x - mean(x)))
      V_T <- V_W + V_B + V_B / m
      SE_pooled <- sqrt(V_T)
      return(SE_pooled)
    }
    ses_pooled <- pool_ses(est, se)
    
    return(data.frame(AME = est_pooled,
                      SE = ses_pooled))
  }
  
  ## Pool estimates
  pooled_imputations <- pool_imputations(list_mfx_summaries)
  
  ## Initialize ourput object
  pooled_summary <- list_mfx_summaries[[1]]
  
  ## Drop unneccesary columns
  pooled_summary$AME <-
    pooled_summary$SE <-
    pooled_summary$z <-
    pooled_summary$p <-
    pooled_summary$lower <-
    pooled_summary$upper <- NULL
  
  ## Pool point estimates and SEs
  pooled_summary$AME <- pooled_imputations$AME
  pooled_summary$SE <- pooled_imputations$SE
  
  ## Add upper and lower bounds
  pooled_summary$lower <-
    pooled_summary$AME + qnorm(.025) * pooled_summary$SE
  pooled_summary$upper <-
    pooled_summary$AME + qnorm(.975) * pooled_summary$SE
  
  
  ## Return
  return(pooled_summary)
}